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ABSTRACT 
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CS| We present a comparison between the observational data on the kinematical 

\ structure of Gl in M31, obtained with the Hubble WFPC2 and STIS instruments, 

^ , and the results of dynamical simulations carried out using the special-purpose 

^ \ computer GRAPE-6. We have obtained good fits for models starting from single 

i \ cluster King-model initial conditions and even better fits when starting our sim- 

^ ! ulations with a dynamically constructed merger product of two star clusters. In 

the latter case, the results from our simulations are in excellent agreement with 
the observed profiles of luminosity, velocity dispersion, rotation, and ellipticity. 
We obtain a mass-to-light ratio of MIL = 4.0 ± 0.4 and a total cluster mass of 
M = (8 ± 1) x 10 6 Mq. Given that our dynamical model can fit all available 

^ | observational data very well, there seems to be no need to invoke the presence of 

an intermediate- mass black hole in the center of Gl. 

Subject headings: black hole physics — globular clusters: individual (Gl) — methods: 

"Eg ! N-body simulations — stellar dynamics 

C3 ■ 



1. Introduction 

We report results from a series of iV-body simulations for the globular cluster Gl in 
M31. Gl is one of the brightest and most massive globular clusters in the local group. Its 
total luminosity (My = —10.94 mag) and central velocity dispersion (ctq = 25. 1 ± 1.7 km/sec) 
are larger than those of any galactic globular cluster (Meylan et al. 2001; Djorgovski et al. 
1997). 
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Meylan et al. (2001) used virial mass estimates and mass estimates from King-Michie 
models for Gl. They obtained total masses in the range 7.3 x 1O 6 M to 15 x 1O 6 M and 
(for their model 4) a core radius, half-mass radius, and tidal radius of 0.53, 13.2, and 187 
pc, respectively. The estimated half-mass relaxation time was 50 Gyr, much longer than the 
Hubble time. 

Gebhardt et al. (2002) have reported evidence for an intermediate-mass black hole of 
2.0lo'g x 1O 4 M in the center of Gl. Based on velocity profiles obtained with the Hubble space 
telescope's STIS instrument, they constructed orbit-based axisymmetric models. Varying 
M/L and the mass of the central black hole, they found a best fit for M/L = 2.5 and 
Mbh = 2 x 10 4 M Q . A model without a central black hole was rejected at a 2a level. 

The presence of such a black hole would be very interesting, for at least two reasons. 
First, it would lie neatly on the Mbh — <* relation for galaxies (Gebhardt et al. 2000; Fer- 
rarese and Merritt 2000). Second, Gl would then be a good example of the type of cluster 
postulated by Ebisuzaki et al. (2001), some of which may find their way into the center of a 
galaxy by dynamical friction, where their intermediate-mass black holes may then merge to 
provide the seeds for supermassive black holes. 

However, before embracing such an exciting conclusion it is all the more important to 
ensure that more conventional explanations of the observational data are ruled out. To this 
end, we have tried to construct the best possible evolutionary model for Gl as a large globular 
cluster that is still in the early stages of core collapse, without harboring an intermediate- 
mass black hole. We have run a set of models with varying initial density profiles, half mass 
radii, total masses, and global M/L until we found a model that gave the best fit to the 
light and velocity profiles of Gl. 

In §2 we describe our numerical method. In §3 we present the results of simulations 
starting with a single non-rotating cluster, and in §4 we show what happens when we consider 
Gl to be the rotating product of a merger of two smaller globular clusters. We briefly 
summarize in §5. 

2. Modeling Method 

In order to model the evolution of Gl using iV-body simulations, we face a scaling and 
a fitting problem: we can only handle ~ 10 5 particles while Gl contains ~ 10 7 stars; and 
we do not know which values to assign to the initial cluster model parameters such as the 
total mass and the half-mass radius. We solve the scaling problem by scaling the dynamical 
parameters in such as way as to reproduce in our model simulations the correct two-body 
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relaxation time scales inferred for Gl from observations. We solve the fitting problem by 
carrying out a large enough number of runs to allow us to isolate simulations that closely 
reproduce the observational data. Without the use of the GRAPE-6 computers (Makino et 
al. 2003), it would have been unpractical to run the several dozen runs needed to determine 
our best fits. 

We used Aarseth's iV-body code NBODY4 (Aarseth 1999). All simulated clusters con- 
tained N = 65,536 stars initially, with a range of masses following Kroupa (2001)'s mass 
function with lower and upper mass limits of 0.1 and 30 M Q , respectively. Our simulations 
did not contain primordial binaries, which is a reasonable simplification for a cluster that 
is still quite far from core collapse. We did not include M31's galactic tidal field, which 
would have a negligible influence at the position of Gl, at least 40 kpc from the center of 
M31. Since tidal effects are unimportant, we are left with two evolution mechanisms: stellar 
evolution and two-body relaxation. 

Stellar evolution was modeled according to the fitting formulae of Hurley et al. (2000), 
using a metallicity of [Fe/H] = —0.95, similar to the mean metallicity of Gl as determined 
by Meylan et al. (2001). We assumed a retention fraction of neutron stars of 15%. 

All simulations were carried out for 13 Gyr and the final density and velocity profiles 
were obtained from 10 snapshots spanning a 500 Myr period centered at T = 12 Gyr. For 
comparison of our models with the observations of Gl, we assume a distance of 770 kpc 
to M31, so one arcsecond corresponds to 3.7 pc. Typically, about 1% of the stars escaped 
from the cluster during a simulation. Only bound stars were used for the comparison with 
observations. 

We have to scale the parameters of our simulations in order to match the most important 
stellar evolution and stellar dynamical parameters of the actual Gl cluster. In order to match 
the relaxation time of Gl, we have to increase the radius of our cluster by 



where subscripts Gl and S denote, respectively, the actual values for Gl and those used 
in our simulations. Effects which strongly depend on the number of particles in the cluster 
are unimportant before the cluster goes into core collapse, so our models should give a valid 
description of the dynamical evolution of Gl up to the present time. 

In the first set of simulations, we started from isotropic King model conditions with no 
initial mass-segregation, and dimensionless central concentrations in the range 4.0 < Wo < 
11.0. For each choice of initial density profile, we ran full simulations for a number of choices 
for the initial physical half-mass radius r h (t = 0) and mass M(t = 0) of Gl until we could 
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fit the surface density profile of Meylan et al. (2001) over a maximum range in radius while 
simultaneously obtaining an optimal fit to the observed velocity profile. For the surface 
density, we scaled our predicted profile by a multiplicative factor to obtain the best fit (in 
practice changing the M/L predicted by our assumed IMF by a factor of 1.5 ~ 2). For the 
velocity profile, we used the symmetrized profile shown by Gebhardt et al. (2002) in their 
Fig. 1 and the ground-based value of Djorgovski et al. (1997), who measured a velocity of 
25.1 ± 1.7 km/sec inside an aperture of 1'.'15 x 7'.'0. For each run, a best fit was determined by 
a x 2 test against the combined data. With improved estimates for r h (t = 0) and M(t = 0), 
a new initial half-mass radius could be calculated and a new simulation was performed. 
Simulations were performed until the half-mass radius changed by less than 5% between 
successive iterations. A more detailed description of our simulations and their results will 
be presented in a forthcoming paper (Baumgardt et al. 2003b). 



3. Single Cluster Simulations 

Figure 1 shows the data for the best fit from among the runs where we started with 
a single cluster in the form of a King model; here the initial central potential depth was 
Wq = 7.5. The top panel shows the inferred projected luminosity density. We can see that 
the fit is very good for r < 15 pc. The reason why the model density drops off sharply at 
large radius is because the initial King model had a tidal radius of 32 pc if we scale it to 
Gl. Two-body relaxation begins to produce an extended halo with a surface density slope 
~ —4, but a Hubble time is too short to let this effect propagate very far into the observed 
halo. Starting from deeper King models (W = 8 or higher) does not solve this problem: 
such models predict too high a surface density around r = 10 pc, while still falling short at 
larger radii. The implication is that Gl must have started with a density distribution more 
extended than any King model that can be fit to the bulk of the stars. 

The bottom panel shows the velocity dispersion inferred from our W = 7.5 model, as 
compared to the dispersion observed by Gebhardt et al. (2002). For larger W values, our 
models produce velocities that are too high at the largest observed radii. Models with slightly 
lower concentration give a somewhat better fit, but when we require the model to reproduce 
the density as well, the combined requirements clearly point to Wq = 7.5 as producing the 
best agreement, and one which falls within the observational errors everywhere except near 
the tidal radius artificially imposed by the initial conditions; we address this limitation in 
the next section. 

Note that our model cluster has a mass smaller than those found from multi-mass King 
model fits by Meylan et al. (2001). Their extreme values stem from the implicit King-model 
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Fig. 1. — The best-fit result, starting from a single W = 7.5 King model. The top panel 
shows the surface luminosity profile S of the simulations (full line) and the observations (filled 
circles). The middle panel shows AS/S Q b s , where AS = S mo d c i — S G b s - The bottom panel 
shows the velocity dispersion in the simulations (solid curve, with dashed curves indicating 
the la error) and observations (filled triangles with error bars). 



- 6- 



requirement that a cluster has complete mass-segregation, which is unphysical in a massive 
cluster like Gl where the relaxation time is much longer than a Hubble time. 

To sum up, an evolutionary model starting from a King model without initial mass 
segregation reproduces both the luminosity profile and the velocity dispersion profile of Gl 
rather well. The fits are not perfect, though, on two counts. First, the best-fit model 
still produces too steep a luminosity profile at larger radii. Second, since we start from a 
spherically symmetric non-rotating model, in principle we cannot fit the observed rotation 
profile or ellipticity. The question is whether we can introduce rotation while simultaneously 
at least preserving, and hopefully improving, the reasonable fits obtained so far. In the next 
section we answer this question affirmatively 

4. Merger Simulations 

Currently favored scenarios for the formation of star clusters are the collapse of giant 
molecular clouds or the collision of smaller clouds (Fall and Rees 1985; Fujimoto and Kumai 
1997). A collision scenario could easily explain the apparent rotation of Gl. It might also 
account for the run of surface density in the halo, since simulations of the merging of two 
stellar systems usually give surface density profiles S(-R) ~ R~ 3 (Sugimoto and Makino 
1989; Makino, Akiyama, & Sugimoto 1990; Okumura et al. 1991). 

Based on these theoretical hints, we have carried out a series of simulations starting with 
an early merger of two star clusters. For the sake of simplicity, we have restricted ourselves 
to the merger of two identical King model clusters on parabolic orbits with initial separation 
r\ = 20 and pericenter separation in A^-body units of p — 1 (another simulation with p = 2 
gave similar results). The chosen pericenter distance corresponds to approximately 1.3 half- 
mass radii for the initial clusters. We used N = 80, 000 equal- mass stars in our merger 
simulation without including stellar evolution, a reasonable approximation given that our 
merger was postulated to occur during formation of the clusters. After the merger product 
had undergone its violent relaxation, we randomly selected 65,536 stars from among all the 
stars still bound to the final merger product, assigned masses drawn from a Kroupa (2001) 
IMF to them, and then started our dynamical evolution simulations for a duration of T = 13 
Gyr. 

Fig. 2 shows the final density and velocity dispersion profiles for our best-fit simulations 
which started with a collision between two Wo = 6.5 initial King models. Note that our 
simulations now reproduce the observed extended halo very well. The agreement between 
the observed and model velocity dispersions is also very good. 
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Fig. 2. — Same as figure 1, but for the merger model that started from two W = 6.5 King 
models. 
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Fig. 3. — Symmetrized radial velocities (top panel) and ellipticity profiles (bottom panel) 
for the observations and our best-fit merger model. The ellipticity is defined as e = f — b/a, 
where a and b are the major and minor axis of the best-fit ellipse for the observations and the 
projected model data. Note that we have omitted the innermost 3 datapoints from Meylan 
et al. (2001), since the ellipticity in the center is not well defined (Meylan 2003, private 
communication) . 
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Figure 3 compares the rotation and ellipticity profiles of our merger model and the ob- 
servations. We measured the rotation profile from two directions perpendicular to each other 
and the minor axis and took the mean of the two directions. The profiles were determined 
from the radial velocities of all bright stars located in an area between angles of 10° and 40° 
with respect to the major axis, in order to make an optimal comparison with Gebhardt et al. 
(2002), who performed their HST/STIS spectroscopy at an angle of 25° against the major 
axis. The agreement between simulated and observed rotation profiles is very good. 

Similarly, the ellipticity profile of our merger run is in good agreement with the ellipticity 
profile of Gl as determined by Meylan et al. (2001), (as can be seen in the lower panel of 
Fig. 3). The iV-body run starts with a near constant ellipticity of about e = 0.25. After 
T = 12 Gyr, the cluster core has become almost spherical due to relaxation effects, while the 
halo ellipticity has remained unchanged. The observations show a similar drop of e toward 
the core. 

Our success in modeling Gl as a merger product does not necessarily imply that a 
merger history is the only way to explain its current state. For example, it is also possible 
that Gl is a heavily stripped remnant of a dwarf spheroidal. What is important is that 
the observed rotation can be well modeled under at least one set of reasonable assumptions, 
as we have shown here. The presence of rotation does not invalidate attempts at modeling 
under the simpler assumption of spherical symmetry; rather it invites a further fine-tuning 
of the already good agreement of spherical models. 

Table 1 summarizes our results for the best-fitting models. It shows W for each initial 
model, as well as the half-mass radius at T = 12 Gyr, and the inferred total mass M and 
M/L required to give the best-fit velocity dispersions. The half-mass radius shown is the 
half-mass radius of our models after they were scaled back to Gl by eq. 1. The errors given 
for the M and M/L values are the statistical errors from the x 2 -^. The global mass-to-light 
ratios we obtain in our best fits lie around M/L ~ 4, relatively large but still within the 
range of mass-to-light ratios observed for galactic globular clusters (Pryor and Meylan 1993). 
Column 6 gives the probability Py that our velocity distribution agrees with the observations 
of Gebhardt et al. (2002) and Djorgovski et al. (1997), determined from a x 2 -test against 
the combined data. 

The last three columns of Table 1 give the Mj L values inside the core (defined as the 
region containing the innermost 1% of bright stars) and the half-mass and core relaxation 
times calculated from the cluster parameters at T = 12 Gyrs and eq. 2-62 of Spitzer (1987). 
Since the half-mass relaxation time is much longer than a Hubble time, Gl has not yet 
reached core collapse and the core of Gl is still dominated by low-mass main-sequence stars. 
Nevertheless, core M/ L values are smaller than global ones since mass- segregation has caused 
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Table 1: Results of the best-fitting iV-body runs. 
Wo Type r h M M/L P v (M/L) c T RH T R (0) 

[pc] [M ] [Gyrs] [Gyrs] 

7.5 Single 6.76 (7.60 ± 0.76) x 10 6 3.80 ± 0.38 0.20 2.41 30.3 0.292 
6.5 Merger 8.21 (8.20 ± 0.85) x 10 6 4.10 ±0.42 0.28 2.44 29.5 0.181 

bright stars, which are more massive than average, to sink into the center. Shortly before 
core collapse, bright stars will be depleted from the center by the heavier neutron stars and 
massive white dwarfs, as in the simulations of M15 of Baumgardt et al. (2003a). 



5. Conclusions 

We have constructed evolutionary models for the massive globular cluster Gl. Starting 
from a Wo = 7.5 King model we can reproduce both the observed luminosity profile for 
r < 15 pc and the observed velocity dispersion profile. A model starting from the merger 
of two Wo = 6.5 King models fares even better: it can reproduce the luminosity, velocity 
dispersion, rotation profiles and ellipticity for the entire range of observations. 

Our simulations were motivated by the recent claim of evidence for a massive central 
black hole. Given that our dynamical model without central black hole can fit all avail- 
able observational data very well, there seems to be no need to invoke the presence of an 
intermediate-mass black hole. Note that we obtained an excellent fit by varying only the 
following basic parameters: the central potential Wq, the initial total mass M(0) and total 
mass-to-light ratio M/L, and the initial half-mass radius r/j(0). Our conclusions are therefore 
robust, and independent of any fine tuning in initial conditions. 

The 2-sigma evidence presented by Gebhardt et al. (2002) for a massive black hole is 
not supported by direct observation of luminosity profile, velocity dispersion and rotation. It 
must have come from the data not presented in their paper (e.g. the higher order moments 
of the velocity profiles, together with multiparameter fits to orbit families). Without inde- 
pendent checks or further observational support, we consider the evidence for such a black 
hole to be inconclusive. 

This work is the first example of the successful detailed dynamical modeling of the 
evolution of a globular cluster with rotation. We have shown how iV-body simulations have 
matured as the most powerful tool to interpret detailed observational data, obviating the 
need for simplifying assumptions such as spherical symmetry or the use of static (e.g. multi- 



- 11 - 



mass King) models. 
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